Indian Subcontinent Rayshader topography map

library(rayshader)

Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

library(raster)
library(tidyverse)

Registered S3 methods overwritten by 'dbplyr':
  method         from
  print.tbl_lazy     
  print.tbl_sql      
-- Attaching packages -------------------------------------------- tidyverse 1.3.1 --
v ggplot2 3.3.4     v purrr   0.3.4
v tibble  3.1.2     v dplyr   1.0.7
v tidyr   1.1.3     v stringr 1.4.0
v readr   1.4.0     v forcats 0.5.1
-- Conflicts ----------------------------------------------- tidyverse_conflicts() --
x tidyr::extract() masks raster::extract()
x dplyr::filter()  masks stats::filter()
x dplyr::lag()     masks stats::lag()
x dplyr::select()  masks raster::select()

library(sf)

Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1

library(stars)

Loading required package: abind

ind_region <- raster::raster("../gt30e060n40.tif")

ind_region_mat <- raster_to_matrix(ind_region)

[1] "Dimensions of matrix are: 4800x6000."

from: https://www.rayshader.com/reference/ambient_shade.html

took around 10 mins to run

ind_region_mat = resize_matrix(ind_region_mat, scale = 2, method = "cubic")

ind_region_mat %>%
 sphere_shade(zscale=3, texture = create_texture("red","red","red","red","white")) %>%
 add_shadow(ambient_shade(ind_region_mat, maxsearch = 100, multicore = TRUE,zscale=1),0) %>%
 add_shadow(lamb_shade(ind_region_mat),0.5) %>%
 plot_map()

Error in { : task 1 failed - "cannot allocate vector of size 879.2 Mb"

Get shape of India

library(rgdal)
library(sp)
library(sf)

ind_shp <- readOGR("V:\\1. R & Python work\\3. R\\India shape file downloaded\\By Survey of India\\India Outline Map\\polymap15m_area.shp")

OGR data source with driver: ESRI Shapefile 
Source: "V:\1. R & Python work\3. R\India shape file downloaded\By Survey of India\India Outline Map\polymap15m_area.shp", layer: "polymap15m_area"
with 314 features
It has 2 fields
Integer64 fields read as strings:  Line_Width 

crop(ind_region, extent(ind_shp))

Error in .local(x, y, ...) : extents do not overlap

myExtent <- spTRansform(ind_shp, CRS(proj4string(ind_region)))

Error in spTRansform(ind_shp, CRS(proj4string(ind_region))) : 
  could not find function "spTRansform"

Tried reducing it to only India from: https://stackoverflow.com/questions/47885065/crop-raster-with-polygon-in-r-error-extent-does-not-overlap

but didn’t work

Trying again from scratch

ind_region <- raster::raster("../gt30e060n40.tif")

ind_region_mat <- raster_to_matrix(ind_region)

[1] "Dimensions of matrix are: 4800x6000."

without resizing this time

ind_region_mat %>%
 sphere_shade(zscale=3, texture = create_texture("red","red","red","red","white")) %>%
 add_shadow(ambient_shade(ind_region_mat, maxsearch = 100, multicore = TRUE,zscale=1),0) %>%
 add_shadow(lamb_shade(ind_region_mat),0.5) %>%
 plot_map()

png("ind_sub_mix_col1.png")
ind_region_mat %>%
 sphere_shade(zscale=3, texture = create_texture("red","red","brown","white","green")) %>%
 add_shadow(ambient_shade(ind_region_mat, maxsearch = 100, multicore = TRUE,zscale=1),0) %>%
 add_shadow(lamb_shade(ind_region_mat),0.5) %>%
 plot_map()
dev.off()

null device 
          1 

# png("ind_sub_mix_col1.png")
ind_region_mat %>%
 sphere_shade(zscale=3, texture = create_texture("brown","red","brown","#77DD77","red")) %>%
 add_shadow(ambient_shade(ind_region_mat, maxsearch = 100, multicore = TRUE,zscale=1),0) %>%
 add_shadow(lamb_shade(ind_region_mat),0.5) %>%
 plot_map()

# dev.off()

# png("ind_sub_mix_col1.png")
ind_region_mat %>%
 sphere_shade(zscale=3, texture = create_texture("#77DD77", "red", "brown", "brown", "white")) %>%
 add_shadow(ambient_shade(ind_region_mat, maxsearch = 100, multicore = TRUE,zscale=1),0) %>%
 add_shadow(lamb_shade(ind_region_mat),0.5) %>%
 plot_map()

# dev.off()

ind_region_mat %>%
 sphere_shade(zscale=3, texture = create_texture("#C4A484", "#C4A484", "#C4A484", "#C4A484","white")) %>%
 add_shadow(ambient_shade(ind_region_mat, maxsearch = 100, multicore = TRUE,zscale=1),0) %>%
 add_shadow(lamb_shade(ind_region_mat),0.5) %>%
 plot_map()

ind_region_mat %>%
 sphere_shade(zscale=3, texture = create_texture("#C4A484", "#C4A484", "#C4A484", "#C4A484","white")) %>%
 add_shadow(ambient_shade(ind_region_mat, maxsearch = 100, multicore = TRUE,zscale=1),0) %>%
 add_shadow(lamb_shade(ind_region_mat),0.5) %>%
 plot_map(rotate = 90)

ind_region_mat %>%
 sphere_shade(zscale=3, texture = create_texture("#C4A484", "#C4A484", "#C4A484", "#C4A484","white")) %>%
 add_shadow(ambient_shade(ind_region_mat, maxsearch = 100, multicore = TRUE,zscale=1),0) %>%
 add_shadow(lamb_shade(ind_region_mat),0.5) %>%
 plot_3d(ind_region_mat)

render_snapshot(filename = "ind_sub_3Dplot7.png")

Trying to have color background but it fails

ind_region_mat %>%
 sphere_shade(zscale=3, texture = create_texture("#C4A484", "#C4A484", "#C4A484", "#C4A484","white")) %>%
 add_shadow(ambient_shade(ind_region_mat, maxsearch = 100, multicore = TRUE,zscale=1),0) %>%
 add_shadow(lamb_shade(ind_region_mat),0.5) %>%
 plot_map(background = "#77DD77")

Error in h(simpleError(msg, call)) : 
  error in evaluating the argument 'x' in selecting a method for function 'plotRGB': unused argument (background = "#77DD77")

Crop with shape

https://www.youtube.com/watch?v=UP2Za1TizOc

ind_outline <- sf::st_read("V:\\1. R & Python work\\3. R\\India shape file downloaded\\By Survey of India\\India Outline Map\\polymap15m_area.shp")

Reading layer `polymap15m_area' from data source 
  `V:\1. R & Python work\3. R\India shape file downloaded\By Survey of India\India Outline Map\polymap15m_area.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 314 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 2815341 ymin: 2177499 xmax: 5678865 ymax: 5444567
Projected CRS: LCC_WGS84

ind_outline %>% 
        st_as_sf() %>% 
        ggplot() +
        geom_sf()

ind_outline

Simple feature collection with 314 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 2815341 ymin: 2177499 xmax: 5678865 ymax: 5444567
Projected CRS: LCC_WGS84
First 10 features:
   Id Line_Width                       geometry
1   0       1875 POLYGON ((5547296 2230982, ...
2   0       1875 POLYGON ((5560180 2232030, ...
3   0       1875 POLYGON ((5549993 2253154, ...
4   0       1875 POLYGON ((5542651 2256150, ...
5   0       1875 POLYGON ((5533962 2260494, ...
6   0       1875 POLYGON ((5523175 2264240, ...
7   0       1875 POLYGON ((3223295 2294948, ...
8   0       1875 POLYGON ((5502051 2315325, ...
9   0       1875 POLYGON ((5522126 2328209, ...
10  0       1875 POLYGON ((5480027 2338995, ...

st_transform(ind_outline, crs = st_crs(ind_outline))

Simple feature collection with 314 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 2815341 ymin: 2177499 xmax: 5678865 ymax: 5444567
Projected CRS: LCC_WGS84
First 10 features:
   Id Line_Width                       geometry
1   0       1875 POLYGON ((5547296 2230982, ...
2   0       1875 POLYGON ((5560180 2232030, ...
3   0       1875 POLYGON ((5549993 2253154, ...
4   0       1875 POLYGON ((5542651 2256150, ...
5   0       1875 POLYGON ((5533962 2260494, ...
6   0       1875 POLYGON ((5523175 2264240, ...
7   0       1875 POLYGON ((3223295 2294948, ...
8   0       1875 POLYGON ((5502051 2315325, ...
9   0       1875 POLYGON ((5522126 2328209, ...
10  0       1875 POLYGON ((5480027 2338995, ...

from: https://r-spatial.github.io/stars/articles/stars1.html

ind_region_stars <- stars::read_stars("../gt30e060n40.tif")

ind_region_stars

stars object with 2 dimensions and 1 attribute
attribute(s), summary of first 1e+05 cells:
                 Min. 1st Qu. Median     Mean 3rd Qu. Max.
gt30e060n40.tif   130     793   1052 1302.186    1648 4795
dimension(s):

plot(ind_region_stars, axes = TRUE)

downsample set to c(10,10)

ggplot() +
        geom_stars(data = ind_region_stars) +
        scale_fill_viridis_c(option = "plasma")

ggplot() +
        geom_stars(data = ind_region_stars) +
        scale_fill_viridis_c() +
        geom_sf(data = ind_outline, alpha = 0)

ind_region_stars_cropped <- st_crop(ind_region_stars, ind_outline) 

Error in st_crop.stars(ind_region_stars, ind_outline) : 
  for cropping, the CRS of both objects have to be identical

st_crs(ind_outline)

Coordinate Reference System:
  User input: LCC_WGS84 
  wkt:
PROJCRS["LCC_WGS84",
    BASEGEOGCRS["WGS 84",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["Degree",0.0174532925199433]]],
    CONVERSION["unnamed",
        METHOD["Lambert Conic Conformal (2SP)",
            ID["EPSG",9802]],
        PARAMETER["Latitude of false origin",24,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",80,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",12.472944,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",35.172806,
            ANGLEUNIT["Degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",4000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",4000000,
            LENGTHUNIT["metre",1],
            ID["EPSG",8827]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]

st_crs(ind_region_stars)

Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    DATUM["World Geodetic System 1984",
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    ID["EPSG",4326]]

from: https://stackoverflow.com/questions/30287065/convert-lambert-conformal-conic-projection-to-wgs84-in-r

library(rgdal)

rgdal: version: 1.5-23, (SVN revision 1121)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 3.2.1, released 2020/12/29
Path to GDAL shared files: C:/Users/vinee/Documents/R/win-library/4.1/rgdal/gdal
GDAL binary built with GEOS: TRUE 
Loaded PROJ runtime: Rel. 7.2.1, January 1st, 2021, [PJ_VERSION: 721]
Path to PROJ shared files: C:/Users/vinee/Documents/R/win-library/4.1/rgdal/proj
PROJ CDN enabled: FALSE
Linking to sp version:1.4-5
To mute warnings of possible GDAL/OSR exportToProj4() degradation,
use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
Overwritten PROJ_LIB was C:/Users/vinee/Documents/R/win-library/4.1/rgdal/proj

crs <- CRS("+proj=lcc +lat_1=30 +lat_2=60 +lat_0=38 +lon_0=126 +datum=WGS84")
ind_outline_crs <- SpatialPoints(ind_outline, proj4string=crs)

Error in (function (classes, fdef, mtable)  : 
  unable to find an inherited method for function ‘coordinates’ for signature ‘"sf"’

from: https://stackoverflow.com/questions/68176438/how-to-adjust-raster-shapefile-projections-in-r-to-make-it-suitable-for-croppi?noredirect=1#comment120495418_68176438

ind_outline <- st_transform(ind_outline, crs = st_crs(ind_region_stars))

ind_outline

Simple feature collection with 314 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 68.14811 ymin: 6.75057 xmax: 97.40683 ymax: 37.08833
Geodetic CRS:  WGS 84
First 10 features:
   Id Line_Width                       geometry
1   0       1875 POLYGON ((93.73132 7.236355...
2   0       1875 POLYGON ((93.84585 7.234451...
3   0       1875 POLYGON ((93.77428 7.429366...
4   0       1875 POLYGON ((93.71204 7.462109...
5   0       1875 POLYGON ((93.63901 7.507875...
6   0       1875 POLYGON ((93.54688 7.550112...
7   0       1875 POLYGON ((73.06285 8.303476...
8   0       1875 POLYGON ((93.40323 8.01913,...
9   0       1875 POLYGON ((93.59233 8.115995...
10  0       1875 POLYGON ((93.22764 8.246983...

ind_outline %>% 
        st_as_sf()

Simple feature collection with 314 features and 2 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 68.14811 ymin: 6.75057 xmax: 97.40683 ymax: 37.08833
Geodetic CRS:  WGS 84
First 10 features:
   Id Line_Width                       geometry
1   0       1875 POLYGON ((93.73132 7.236355...
2   0       1875 POLYGON ((93.84585 7.234451...
3   0       1875 POLYGON ((93.77428 7.429366...
4   0       1875 POLYGON ((93.71204 7.462109...
5   0       1875 POLYGON ((93.63901 7.507875...
6   0       1875 POLYGON ((93.54688 7.550112...
7   0       1875 POLYGON ((73.06285 8.303476...
8   0       1875 POLYGON ((93.40323 8.01913,...
9   0       1875 POLYGON ((93.59233 8.115995...
10  0       1875 POLYGON ((93.22764 8.246983...

ggplot() +
        geom_stars(data = ind_region_stars) +
        scale_fill_viridis_c() +
        geom_sf(data = ind_outline, alpha = 0)

ind_region_stars_cropped <- sf::st_crop(ind_region_stars, ind_outline) 

Error in s2_geography_from_wkb(x, oriented = oriented, check = check) : 
  Evaluation error: Found 1 feature with invalid spherical geometry.
[1] Loop 67 is not valid: Edge 31 is degenerate (duplicate vertex).

box = c(xmin = 68.14811, ymin = 6.75057, xmax = 97.40683, ymax = 37.08833)

plot(st_crop(ind_region_stars, box))

Error in h(simpleError(msg, call)) : 
  error in evaluating the argument 'x' in selecting a method for function 'plot': no applicable method for 'st_geometry' applied to an object of class "c('double', 'numeric')"

ind_region_stars_cropped <- sf::st_crop(ind_region_stars, ind_outline, crop = FALSE) 

Error in s2_geography_from_wkb(x, oriented = oriented, check = check) : 
  Evaluation error: Found 1 feature with invalid spherical geometry.
[1] Loop 67 is not valid: Edge 31 is degenerate (duplicate vertex).

from: https://www.youtube.com/watch?v=xbtyaja8tro&t=225s

plot(ind_region)
plot(ind_outline, add = TRUE)

Warning in plot.sf(ind_outline, add = TRUE) :
  ignoring all but the first attribute

ind_region_raster_cropped <- raster::crop(ind_region, ind_outline)
ind_region_raster_cropped

class      : RasterLayer 
dimensions : 3641, 3511, 12783551  (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : 68.15, 97.40833, 6.75, 37.09167  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : gt30e060n40 
values     : 1, 8752  (min, max)

ind_region_raster_masked <- mask(ind_region, ind_outline)
ind_region_raster_masked

class      : RasterLayer 
dimensions : 6000, 4800, 28800000  (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : 60, 100, -10, 40  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : gt30e060n40 
values     : 1, 8238  (min, max)

ind_region_raster_crop_masked <- mask(ind_region_raster_cropped, ind_outline)
ind_region_raster_crop_masked

class      : RasterLayer 
dimensions : 3641, 3511, 12783551  (nrow, ncol, ncell)
resolution : 0.008333333, 0.008333333  (x, y)
extent     : 68.15, 97.40833, 6.75, 37.09167  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : memory
names      : gt30e060n40 
values     : 1, 8238  (min, max)

write tif

writeRaster(ind_region_raster_crop_masked, filename="India_region.tif", format="GTiff", overwrite=TRUE)

plotting India with rayshader

ind_only_mat <- raster_to_matrix(ind_region_raster_crop_masked)

[1] "Dimensions of matrix are: 3511x3641."

ind_only_mat %>%
 sphere_shade(zscale=3, texture = create_texture("#C4A484", "#C4A484", "#C4A484", "#C4A484","white")) %>%
 add_shadow(ambient_shade(ind_only_mat, maxsearch = 100, multicore = TRUE,zscale=1),0) %>%
 add_shadow(lamb_shade(ind_only_mat),0.5) %>%
 plot_map()

colors used from: https://www.google.com/search?q=india+relief+map&client=firefox-b-d&sxsrf=ALeKk01qdHNINMtsDP-Yr1uL_bcvzuphoA:1624969260878&tbm=isch&source=iu&ictx=1&fir=eHLzW_TeaFkWrM%252C8wqDRm6USS1W9M%252C_&vet=1&usg=AI4_-kRg3bVBakksvo7YtWKASkgUWV7AEw&sa=X&ved=2ahUKEwiI0MKc6rzxAhUVQH0KHZMPA2kQ9QF6BAgEEAE#imgrc=eHLzW_TeaFkWrM

ind_only_mat %>%
 sphere_shade(zscale=3, texture = create_texture("#C4A484", "#C4A484", "#C4A484", "#C4A484","#e7c6b2")) %>%
 add_shadow(ambient_shade(ind_only_mat, maxsearch = 100, multicore = TRUE,zscale=1),0) %>%
 add_shadow(lamb_shade(ind_only_mat),0.5) %>%
 plot_map()

colors used from: https://www.google.com/search?q=india+relief+map&client=firefox-b-d&sxsrf=ALeKk01qdHNINMtsDP-Yr1uL_bcvzuphoA:1624969260878&tbm=isch&source=iu&ictx=1&fir=dpw46_M3IQqdSM%252CV25s7No3fn6eyM%252C_&vet=1&usg=AI4_-kTVqgKmNiQE8qgYQTiuK52KyrioPA&sa=X&ved=2ahUKEwiI0MKc6rzxAhUVQH0KHZMPA2kQ9QF6BAgHEAE#imgrc=dpw46_M3IQqdSM

ind_only_mat %>%
 sphere_shade(zscale=3, texture = create_texture("#e7c6b2", "#e7c6b2", "#e7c6b2", "#e7c6b2","#b0b685")) %>%
 add_shadow(ambient_shade(ind_only_mat, maxsearch = 100, multicore = TRUE,zscale=1),0) %>%
 add_shadow(lamb_shade(ind_only_mat),0.5) %>%
 plot_map()

ind_map_green <- ind_only_mat %>%
 sphere_shade(zscale=3, texture = create_texture("#e7c6b2", "#e7c6b2", "#e7c6b2", "#e7c6b2","#e0f0e3")) %>%
 add_shadow(ambient_shade(ind_only_mat, maxsearch = 100, multicore = TRUE,zscale=1),0) %>%
 add_shadow(lamb_shade(ind_only_mat),0.5) %>%
 plot_map()

typeof(ind_map_green)

[1] "NULL"

ind_only_mat %>%
 sphere_shade(zscale=3, texture = create_texture("#e7c6b2", "#e7c6b2", "#e7c6b2", "#e7c6b2","#779ecb")) %>%
 add_shadow(ambient_shade(ind_only_mat, maxsearch = 100, multicore = TRUE,zscale=1),0) %>%
 add_shadow(lamb_shade(ind_only_mat),0.5) %>%
 plot_map()

ggplot() +
        labs(title = "India") +
        theme(plot.background = element_blank(),
              panel.background = element_blank())

library(magick)

Linking to ImageMagick 6.9.12.3
Enabled features: cairo, freetype, fftw, ghostscript, heic, lcms, pango, raw, rsvg, webp
Disabled features: fontconfig, x11

image_with_text <- image_annotate(base_img, "India", size = 70, color = "#779ecb",
                                  location = "+450+120") %>% 
        image_annotate(., "created by: ViSa", color = "grey",size = 15, location = "+490+730")

image_write(image_with_text, "ind_brown_withtext_blue.png")